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We use the kinetic theory of gases to compute the Kolmogorov-Sinai en¬ 
tropy per particle for a dilute gas in equilibrium. For an equilibrium system, 
the KS entropy, hxs is the sum of all of the positive Lyapunov exponents 
characterizing the chaotic behavior of the gas. We compute where N 

is the number of particles in the gas. This quantity has a density expansion of 
the form h^s/AI = on[— lnh+6+0(h)], where n is the single-particle collision 
frequency and h is the reduced number density of the gas. The theoretical 
values for the coefficients a and b are compared with the results of computer 
simulations, with excellent agreement for a, and less than satisfactory agree¬ 
ment for b. Possible reasons for this difference in b are discussed. 

One of the important quantities characterizing the chaotic behavior of a dynamical sys¬ 
tem is the Kolmogorov-Sinai (KS) entropy, h-KS [0)3- If fhe system is isolated and closed, i.e. 
there is no escape of particles from the system, then according to Pesin’s theorem 0], h^s 
is the sum of all the positive Lyapunov exponents of the system, where the Lyapunov ex¬ 
ponents characterize the rate of exponential separation of the system’s trajectories in phase 
space. This quantity is of considerable interest because a positive KS entropy implies that 
the system is chaotic, and because the KS entropy measures the rate at which information 
about the phase space trajectories is produced by the dynamics. Further, sums of Lyapunov 
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exponents figure prominently in the expressions for transport coefficients of fluids in terms 
of dynamical quantities Recently methods based upon the kinetic theory of gases 

have been applied to compute the chaotic properties of simple systems such as hard disk 
and hard sphere Lorentz gases [^-|^, as well as Lorentz lattice gases This work has 

shown that it is possible to calculate theoretically Lyapunov exponents for Lorentz gases at 
low density under a variety of equilibrium and non-equilibrium situations. The theoretical 
results agree very well with computer simulations for the same quantities. Here we show 
how kinetic theory methods can be extended to a calculation of the equilibrium KS entropy 
of dilute gases composed of N particles interacting with short range forces. This extension 
of the methods developed earlier for Lorentz gases to systems of N moving particles allows 
us to begin a thorough exploration of the chaotic properties of systems with large numbers 
of moving particles, as well as of Lorentz gases, both in and outside equilibrium states. In 
subsequent publications Van Beijeren and Van Zon will describe methods to obtain the 
largest positive Lyapunov exponent for a dilute gas in equilibrium. 

We consider a gas composed of N identical particles in equilibrium. The particles obey 
classical mechanics, each have mass m = 1, and their positions and velocities are denoted 
by Tj, Uj, respectively, with i = 1, The particles interact with strongly repulsive, short 

range forces with a hnite range of the interparticle force, a. In the case that the particles 
are hard disks or hard spheres, a is the diameter of each particle. Our goal is to determine 
the KS entropy per particle, kxs/^, iii the thermodynamic limit, assuming that the system 
is in thermal equilibrium at temperature T, and at low densities with na'^ -C 1, where n is 
the number density, n = N/V, V is the volume of the system, and d is the dimension of the 
system, d = 2,3. 

We consider the trajectory of our system on the constant energy surface in 
the 2(iiV-dimensional phase space, T, starting from some initial point Xq = 
(ri(0), ni(0), r2(0), h2(0),..., rAr(O), h7v(0)). We then consider a bundle of trajectories that 
start at inhnitesimally nearby points in phase space, with each trajectory in the bundle 
denoted by X(0) -|-5X(0), for some inhnitesimal 5X(0) = (5ri(0), 5ni(0),..., 5rAr(0), ^^^(0)). 
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The time rate of separation of this trajectory bundle in the various directions perpendicu¬ 
lar to the direction of flow, if exponential, is characterized by a set of non-zero Lyapunov 
exponents which are positive in the expanding directions and negative in the contracting 
ones. The KS entropy is, for our system, given by the sum of all of the positive Lyapunov 
exponents, by Pesin’s theorem [|^^. In order to determine the Lyapunov exponents and 
their sum, we need to follow the dynamics of the trajectory bundle over a time t very large 
compared to some characteristic microscopic time, to, which for low density systems can be 
taken to be the mean free time between collisions for a typical particle. Recently Dellago, 
Posch and Hoover, as well as Gaspard and Dorfman have obtained the equations of motion 
for a trajectory bundle in phase space for a system of hard disk or hard sphere particles 
These equations also apply to any system of particles with strong short range 
interaction potentials, provided we neglect effects due to collisions involving more than two 
particles at a time, or due to the hnite duration of collisions taking place in the fluid. To 
lowest order in the density, these multiple collision and hnite collision time effects may be 
ignored and the dynamics may be thought of as periods of free motion of the particles punc¬ 
tuated by binary collisions between pairs of particles. During the free motion, the positions 
and velocities of the particles change with time according to: 


Tiit) = ri(0) + Vi{0)t 

(1) 

Vi{t) = Ui(0) for i = 

(2) 


Here the time t = 0 represents the instant of the last binary collision in the system, and 
time t is some time between the last binary collision and the next. When a collision takes 
place between particles k and I, say, there is an effectively instantaneous change in the 
velocities of particles k, I while the positions and velocities of all of the other N — 2 particles 
are unaffected. Immediately after the {k, 1) collision, the positions and velocities of the two 
particles, k, I are given by 

4 = + [vik ■ n]h, 

v'i = vi- [vik ■ n\h, 
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( 3 ) 

( 4 ) 





fi = fk + ah. 


( 5 ) 


Here n is a unit vector from the center of particle k to the center of particle I at the point of 
closest approach during the binary collision, vik = vi — Vk, and vik ■ h < 0. The free motion 
of the other N — 2 particles i ^ k,l is not affected by the k, I collision. Now we consider the 
other trajectories in the nearby bundle. We assume that the bundle is sufficiently narrow 
that all trajectories in the bundle exhibit the same collisions in the same sequence, with 
slight differences in the times of the collisions, in the velocities and positions before and 
after the collisions, and in the points of closest approach. The analysis of Dellago, Posch, 
and Hoover leads immediately to equations of motion for the deviations of the positions and 
velocities in the trajectory bundle from the main trajectory. During the free motion between 
collisions, (jrj, and 6vi evolve according to 


6ri(t) = (jfi(O) + t6vi{0) 

(6) 

Svi(t) = Svi{0). 

(7) 


Whenever a collision takes place in the system, the deviations in positions and velocities of 
the colliding particles are directly affected, while those of the other particles are not. The 
velocity deviations for the colliding pair, say particles k, I, can be obtained from Eqs. (3,4) 
by taking linear deviations, that is, 

= + ( 8 ) 
5v'ik = 5v'i- 5v'k = Svik - 2{Svik ■ h)h - 2[{vik ■ Sh)h + {vik ■ h)Sh], (9) 

where for later convenience we have introduced the center of mass velocity Vik = (1/2) [f4 + 
vi], in addition to the relative velocity of the colliding pair of particles. Here Sh is the 
inhnitesimal displacement of the unit vector in the direction of closest approach, due to the 
displacement of the trajectory, and h- 6h = 0. To calculate 6h we use the coupled equations 
for the location of the collision for both the nndisplaced and the displaced trajectories: 

o-h = fi^o + nvi - fkfi - TkVk 

a6h = Srifi - 6rkfi + vikSnk + nSvi - TkSvk, 
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( 10 ) 

( 11 ) 



where we have kept terms to linear order in the deviations. Here r;^o is the position of particle 
I at the instant of its last collision with another particle in the system, and similarly for r^. q- 
Also Ti is the time elapsed from the last collision of particle I until its collision with k, with 
a similar dehnition for Finally Srik is the difference in time between the I, k collision for 
the displaced and the undisplaced trajectories. We assume here that in our calculation of 
hxs/N, the dominant terms will come from terms proportional to Tki = {rk + r/)/2 which is 
on the order of the free time between collisions of particles in the system and which scales 
inversely with the density of the gas. For low densities this free time can be quite large, and 
our assumption is that we can neglect the term 5rifl — 5rkfl in Eq. (|^) compared to terms 
proportional to the times between collisions for each of the particles. This assumption will 
be discussed further in the concluding remarks. Then using the condition h • (5h = 0, and 
changing to center-of-mass and relative velocities, we hnally obtain equations for the change 
in the displaced center of mass and displaced relative velocities of particles k, I at the k, I 
collision: 


6v[ = 6vi for i ^ k,l 

= SVki 

5v'ik = 6vik - 2{6vik ■ n)n 


( 12 ) 

(13) 




{vik ■ 6Vki) - ^A 5Vki ■n)]h + 5Vki{vik ■ n) - vik{SVki ■ h) 

[vik ■ n) 


T, 


kl 


Vik ■ Svik) 


■'Ik 


{vik ■ n) 


{6vik • h) h + {vik ■ n)Svik - vik{5vik ■ fi) 


(14) 


Here tik = Ti — Tk- Using these equations we can express the change of the velocity deviations 
of all of the particles at a k,l collision as a matrix M^,; acting on the column vector 6v whose 
elements are the velocity deviations, 6vi for ^ = 1, as ■ 6v. Here has 

the value 1 along the diagonal for the elements corresponding to the N — 2 particles which 
are not involved in the k, I collision, and non-zero elements for the changes in the velocity 
deviations for particles k and /, which can be obtained from Eqs. (13, 14) in terms of the 
pre-collision velocity deviations and the times Tki,tik- All other matrix elements of are 


5 








zero. 


Thus by neglecting all of the deviations in positions at the instant of binary collisions 
Sri^o, we can easily write an evolution equation for the velocity deviation vector as the result 
of a sequence of binary collisions in the system as 


Sv(t) = ■ M „2 ■ ■ ■ • 5v(0). 


(15) 


In Eq. (|l^ we have used the facts that the dynamics of the system is a sequence of binary 
collisions, labeled by the subscripts aj with aK(t) the last collision up to time t, where 
only the velocities of two particles change, separated by free particle motions where none 
of the velocities change. For low densities, the product of matrices on the right hand side 


of Eq. (pISD can be thought of as a product of random dN x dN matrices, since there 
are no correlations between the collisions in the sequences. The random elements of the 
sequence are the particles involved in the individual collisions, the collision parameters of 
each collision, and the time intervals, for each particle, between the collisions that it suffers. 
Under these circumstances all of the Lyapunov exponents and the KS entropy of the system 
can be obtained, in principle, by determining the eigenvalues of this product of random 
matrices, using the known distributions of free times and collision parameters for a dilute 
gas in equilibrium. The positive Lyapunov exponents and the KS entropy can be obtained 
by using the fact that almost all trajectories in phase space will separate with time t, as t 
grows large, since the probability of finding two nearby, but otherwise randomly selected, 
trajectories that approach each other for arbitrarily long times, is vanishingly small. The 
negative Lyapunov exponents are obtained by considering the time reversed motion and 
using the fact that trajectories that approach each other in the forward time direction will 
separate in the time reversed motion. This analysis was applied to the random, dilute, three 
dimensional Lorentz gas by Latz, van Beijeren, and Dorfman who were able to calculate 
the positive and negative Lyapunov exponents by analyzing the products of random matrices 
similar to those in Eq. ( 0 ) along these lines. Due to the increased size of the matrices 
considered here, the determination of the individual Lyapunov exponents is still a formidable 
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analytical problem. However the determination of the KS entropy is relatively elementary 
as we now show. 

Snppose that the eigenvalnes of the prodnct of the matrices on the right hand side of 
Eq. (^) have the form exp(fAi). Then the logarithm of the determinant of this matrix 
prodnct will have the form f ^ A*. Since almost all trajectories will lead to separation, only 
the positive exponents will appear in this snm. Fnrther, if we assnme that the system is 
ergodic then for long times t, all possible collision parameters and free times will appear in 
the matrices. Tims we can write 

^ j=K(t) 

hxs = lim - In I det M„ | 

t^oc t ^ ^ 

Nv 

= — (In I detMi, 2 |). (16) 


where v is the average collision freqnency per particle, so that the nnmber of collisions 
taking place in the gas per nnit time is Ni'/2 (each collision involves two particles), and 
the angnlar brackets denote averages over the collision rates and parameters, the free time 
distribntions, and the velocity distribntions for a dilnte gas in eqnilibrinm. We have also 
nsed the indices 1, 2 to label the colliding particles in a typical matrix M. This matrix has a 
particnlarly simple strnctnre when expressed in terms of the deviations of the center-of-mass 
and relative velocities of particles 1, 2 and the velocity deviations of particles 3,..., iV which 
are not affected by the 1, 2 collision. 


1 0 0 ^ 


Ml,2 = ABO (17) 

\ ° ^ / 

Here we have organized the velocity deviation vector so that = {SVu, Sv 2 i, ..., 5h)v)^, 
where T denotes a transpose. The fonr npper-left snbmatrices are all d x d dimensional. 
The d X d dimensional snbmatrices A, B, are easily obtained from Eqs. (13, 14) above. One 
easily sees that det Mi ,2 = det B, which is simple to calcnlate. Therefore we hnd that 


u 


hxs/N = — (In [ 2 Ti 2 |n 2 i|/(c’'cosd)]) for d = 2 and 


(18) 
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(19) 


= z/ (In I [ 2 Ti 2 \v 2 i\/c] I) for d = 3. 

Here 9 is the angle of incidence in a binary collision, with —7r/2 <9 <7^/2 for two dimensions 
and 0 < 0 < 7r/2 in three dimensions. Eqs. (18, 19) are the central theoretical results of this 
work. We have obtained explicit formulae for the KS entropy for a dilute gas of particles 
with strong short range forces. We now give the explicit form of the integrals required for 
the calculation of hxs for d = 2,3. For d = 2 we obtain 

^ 2^ J J ^^20o(fb)|h'i2| J d9a{TT - 29)la 

/*00 POO 

xiy{vi)iy{v 2 ) dTi j In [ 2 Ti 2 |'hi 2 |/(acos 0 )]. (20) 

Jo Jo 

Here 4>o{v) is the Maxwell-Boltzmann velocity distribution function, u{v) is the collision 
frequency for particles with velocity v, a{9) is the differential cross section for scattering 
under an angle 9 and J 2 is a normalization factor obtained by replacing the logarithmic term 
in the numerator by unity and equating the resulting expression to unity. The corresponding 
expression for three dimensions is 

j dFi0o(Fi) j ^^ 200 (^ 2 ) 1 ^ 12 ! j d9 sm 9 cos 9a{n - 29) j dcf) 

POO POO 

xiy{v,)iy{iy{v 2 ) / dr, / ln[ 2 Ti 2 |Fi 2 |/a] (21) 

Jo Jo 

To evaluate these expressions we need the velocity dependent collision frequencies for 
particles with velocities Vi,V 2 , since these frequencies are needed for the evaluation of the 
average of lnTi 2 . For particles with interactions of hnite range expressions for these collision 
frequencies are easily obtained from elementary kinetic theory, and explicit results can be 
given for hxs, after some straightforward numerical integrations. We find that for hard 
disks {d = 2) 

Hks/N = ^ [- In(na^) + 0.209 + 0{n)] (22) 

where z/ = [(27r^/^n(T)/(/5m)^/^] is the average collision frequency of a dilute gas of hard disks 
at temperature T = (/cs/?)"^ , and ks is Boltzmann’s constant. For a gas of hard spheres 
we hnd 
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hKs/N = i^[— ln(7m(j^) + 0.562 + 0{n)] 


(23) 


where v = [(47r^/^n(T^)/(/3m)^/^]. 

To check these results we have carried out extensive numerical simulations to compute 
the whole set of Lyapunov exponents for equilibrium hard-disk and hard-sphere systems. 
hxs is obtained as the sum of all the positive exponents. In addition to the phase-space 
trajectory, the “exact” molecular-dynamics (MD) algorithm 0.0 follows the time evolu¬ 
tion of a complete set of inhnitesimal vectors (5X;, / = 1, • • •, 2dN, in tangent space by using 
Eqs. (6) and (7) for the streaming between collisions, and exact collision maps analogous to 
Eq. i) for the discontinuous jumps at a collision. The tangent vectors are periodically re- 
orthonormalized, and the Lyapunov exponents are obtained as the time-averaged logarithms 
of the renormalization factors. 

The MD-method requires the precise localization of the particle collision points and be¬ 
comes inefficient for very low densities. We have therefore applied also a variant of the direct- 
simulation Monte Carlo (DSMC) technique to the computation of the Lyapunov exponents 
I^, where the phase-space dynamics of the system is modeled by appropriate probability 


distributions of the collision parameters and of the collision time, and the tangent-space 
dynamics remains completely deterministic [Q. For very low densities the DSMC-algorithm 
becomes exact and is much more efficient than the MD method. It is expected to fail for 
high densities for which correlated collisions become important. 

The simulations, performed with both methods, are for 36 and 64 disks in two dimensions, 
and for 32 and 108 spheres in three. The usual periodic boundary conditions are employed. 
To present our data we use reduced units for which the particle diameter cr, the kinetic energy 
per particle, K/N = m^-vf/2N, and the particle mass m are unity. Time is measured in 
units of (ma'^N/, and the number density n = N/V in units of where V is the 
volume of the simulation box. n is varied between 10“® and 0.1. Our choice for the unit of 
energy corresponds to a temperature T for which ksT = K/N = 1 in two dimensions, and 
ksT = {2/3)K/N = 2/3 in three. The dependence of kxs/iNu) on the collision frequency 
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V is displayed in Figs. |l| and for the respective two- and three-dimensional systems studied 
here. As indicated by the labels, the diamonds and squares denote MD results, whereas the 
plus signs and crosses refer to DSMC data. It is obvious from the Figures, that there is very 
good agreement between the MD and DSMC results for kxs/Nu, and that this agreement 
persists even to the largest collision frequencies of interest. We therefore do not distinguish 
in the following discussion between these two independent sets of simulation data. 

The particle-number dependence of the Kolmogorov-Sinai entropy per particle is found 
to be very small for hard disks, with the points for 36 and 64 particles practically undistin- 
guishable on the scale of Fig. |I[ For hard spheres there is a noticeable difference between 
the results for N = 32 (dashed) and N = 108 (dotted) in Fig. ||. However, it is expected 
that the larger of the two systems is already close to the thermodynamic limit, and that a 
comparison of our hnite-particle data with the theoretical prediction, valid for N oo, is 
meaningful. 

To facilitate this comparison, we rewrite Eqs. (p^ and (p3|) according to 


hKs/{Nu) = a[-ln(i//i/o) + b] + 0(z/^), (24) 

where Uq = 2{'n /for hard disks, and z/q = 4(7rmcT^)“^/^(2iC/3A^)^/^ for 
hard spheres. Neglecting the higher-order terms and treating a, b as £t parameters, this 
expression is htted to the experimental data for kxs/i^^) range 10“’^ < v < 10“^. 

The results for a and b are summarized in Table |. The theoretical expectation for h^s/iNv) 
is indicated by the smooth line in Fig. |l|, {a,b} = {0.5,0.209}, for disks, and in Fig. |^, 
{a,b} = {1,0.562}, for spheres. 

From an inspection of Table | we conclude that the theoretically expected value for a, 
which determines the slope of the straight lines in the hgures, is almost perfectly reproduced 
by the hts. The theoretical values for the other parameter b, however, are not recovered by 
the htting procedure. This is reflected in the offset of the theoretically expected smooth line 
with respect to the experimental results in the Figs. and |^. 

To shed some light into the origin of this disagreement we have numerically evaluated 
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Eq. (p!6t) by averaging over a molecular-dynamics simulation run. We find perfect agreement 
of this average with the theoretical predictions of Eqs. (^) and (pdf). Moreover, a computer 
simulation where the terms 6fkfi and 6fk,o are set equal to zero for each k, I collision yields 
results that do agree with Eqs. (22, 23), but not with the results where these spatial 
displacements are taken into account. Thus, there appears to be a problem with the neglect 
of the spatial deviations Sri^ — 5rkfl in going from Eq. (11) to Eq. (14). That is, the neglect 
of 5rifi — 5rkfi in Eq. (11) has not yet been justified by an analytic calculation where such 
terms are included. In fact, a more detailed analysis reveals that indeed the inclusion of the 
neglected terms yields corrections to the coefficient b, which may crudely be estimated to 
be close to In 2. However these are very preliminary results which must be verified by more 
extensive calculations. It is worth noting that in the case of the two and three dimensional 
Lorentz gas the neglect of spatial deviations in the equivalent of Eq. |T^ can be fully justified 
and there the theoretical results for the analogs of the coefficients a and b agree very well 
with the simulation results P-|8|,p!^. Finally notice that the assumptions made do not affect 
the theoretical values of the leading term involving only the parameter a. Presently we 
are working both on calculations of the terms neglected above and on a more systematic 
approach using methods based upon the BBGKY hierarchy equations, which will provide a 
systematic density expansion of kxs- 

In conclusion, we have obtained the first theoretical values for the leading term in the 
KS entropy per particle for dilute gases with short range forces, both for two and three 
dimensional systems. The comparison of these predictions with results of computer simu¬ 
lations are excellent, but we do not yet have a final analytic expression for the first order 
corrections to these results. The computer simulations, however, do provide accurate values 
for these corrections, to which eventually analytic results can be compared. There seems to 
be little doubt that our results can be extended to dense gases and to non-equilibrium states 
of interest for the escape-rate or the Gaussian thermostat approaches for relating transport 
coefficients to the chaotic properties of fiuid systems. 
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TABLE I. Fit parameters for the fit of Eq. (23) to the experimental hxs /for a system 
containing N hard disks, (d = 2), or spheres, (d = 3). The fitting range is 10“^ < v < 10“^, where 
u is the collision frequency in the reduced units introduced in the main text. MD and DSMC refer 
to the numerical method used for the computation. The quoted errors for the fit parameters a and 
b are standard deviations. 


d 

N 

Method 

a 

b 

2 

36 

MD 

0.499 ± 0.001 

1.339 ± 0.003 

2 

36 

DSMC 

0.500 ± 0.001 

1.326 ± 0.006 

2 

64 

MD 

0.499 ± 0.001 

1.366 ± 0.005 

2 

64 

DSMC 

0.500 ± 0.001 

1.358 ± 0.009 

2 

oo 

Theory 

0.5 

0.209 

3 

32 

MD 

0.998 ± 0.001 

1.42 ± 0.01 

3 

32 

DSMC 

1.000 ± 0.001 

1.38 ± 0.01 

3 

108 

MD 

1.03 ± 0.03 

1.34 ± 0.29 

3 

108 

DSMC 

1.036 ± 0.002 

1.34 ± 0.02 

3 

oo 

Theory 

1.0 

0.562 
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FIGURE CAPTIONS 


1- hxs/i^^), ^ function of the collision frequency, for planar systems containing 36 

and 64 hard disks, u is measured in the reduced units introduced in the main text. 
The points are labeled according to the computational method, MD or DSMC, the 
dimension, d = 2, and the number of particles N. The smooth line is the theoretical 
prediction according to Eq. and the (almost undistinguishable) dashed and 

dotted lines are hts of Eq. (P^) to the experimental points for iV = 36 and N = 64, 
respectively. The htting range is 10“^ < v < 10“^. The £t parameters a and b are 
listed in Table |I[ 

2. bxs/iNu), as a function of the collision frequency, for three-dimensional systems con¬ 
taining 32 and 108 hard spheres, u is measured in the reduced units introduced in 
the main text. The points are labeled according to the computational method, MD or 
DSMC, the dimension, d = 3, and the number of particles N. The smooth line is the 
theoretical prediction according to Eq. (^), and the dashed and dotted lines are fits 
of Eq. (p^) to the experimental points for N = 32 and N = 108, respectively. The 
htting range is 10“^ < v < 10“^. The £t parameters a and b are listed in Table 
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Fig. 1: van Beijeren et. al. 





Fig. 2: van Beijeren et. al. 



